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Abstract 


The planetary exospheres are poorly known in their outer parts, since the neutral densities 
are low compared with the instruments detection capabilities. The exospheric models are thus 
often the main source of information at such high altitudes. We present a new way to take into 
account analytically the additional effect of the radiation pressure on planetary exospheres. In a 
series of papers, we present with an Hamiltonian approach the effect of the radiation pressure on 
dynamical traje ctories, dens it y profiles an d escap ing thermal flux. Our work is a generalization 


of the study by 


Bishop and Chamberlainl (119891) . In this first paper, we present the complete 


exact solutions of particles trajectories, which are not conics, un der the influence of th e solar 
radiation pressure. This problem w a s rece ntly partly solved by Lantoine and Russelll ( 2011 1 
and completely by Biscan i and Izzo ( 201 4). We give here the full set of solutions, including 
solutions not previously derived, as well as simpler formulations for previously known cases and 


comparis ons wi t 
problem (IStarkl. 


i rece nt works. The solutions given may also be applied to the classical Stark 


191411 : we thus provide here for the first time the complete set of solutions for 


this well-known effect in term of Jacobi elliptic functions. 
Keywords: exosphere, radiation pressure, Stark effect, trajectories 


1. Introduction 


The exosphere is the upper layer of any planetary atmosphere: it is a quasi-collisionless 
medium where the particle trajectories are more dominated by gravity th an by c ollisions. 
Abov e the exobase, the lower limit of the exosphere, the Knudsen number (iFerziger and Kaperl . 
19721 ) becomes large, collisions become scarce, the distribution function cannot be considered 
as maxwellian anymore and, gradually, the trajectories of particles are essentially determined 
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by the gravitation and radiation pressure by the Sun. The trajectories of particles, subject to 


the gravitational force, are comp 


case with the radiation pressure ([Bishop and Chamberla inl. [19891 ) 


ete l y solved w ith the equ ation s of motion, but it is not the 


In the absence of radiation pressure, we can distinguish three types of trajectories for the 
exospheric particles : 

- the escaping particles come from the exobase and have a positive mechanical energy: they 
can escape from the gravitational influence of the planet with a velocity la rger t han the 


1910). They 


escape velocity. These particles are responsible for the Jeans’ escape fjj eansl . 
can also be defined as crossing only once the exobase, 

the ballistic particles also come from the exobase but have a negative mechanical energy, 
they are gravitationally bound to the planet. They reach a maximum altitude and fall 
down on the exobase if they do not undergo collisions. They cross the exobase twice, 
the satellite particles never cross the exobase. They also have a negative mechanical 
energy but their periapsis is above the exobase: they orbit along an entire ellipse around 
the planet without crossing the exobase. The satellite particles result in their major part 


from ballistic particles imdergoing few collisions mainly near the exobase flBeth et ah 
(2014))- Thus, they do not exist in a collisionless model of the exosphere. 

The radiation pressure disturbs the conics (ellipses or hyperbolas) described by the parti¬ 


cles under the influence of gravity. The resonant scattering o 


momentum transfer from the photon to the atom or molecule flBurns et al.l . [l979|). In the non 


solar photo ns leads to a total 


relativistic case, assuming an isotropic reemission of the solar photon, this one is absorbed in 
the Sun direction and scattered with the same probability in all directions. For a sufficient flux 
of photons in the absorption wavelength range, the reemission in average does not induce any 
momentum transfer from the atom/molecule to the photon. The momentum variation, each 


secon d, between bef ore and a : 


’ter t 


re scattering imparts a force, the radiation pressure. 


Bishop and Chamberlain (119891) proposed to analyze its effect on the structure of planetary 


exospheres. In particular, they highlighted analytically the “tail” phenomenon at Earth: the 
density for atomic Hydrogen, which is sensitive to the Lyman-a photons, is higher in the night- 
side direction than in the dayside direction in the Earth corona. Nevertheless, their work was 
limited only to the Sun-planet axis, with a null component assumed for the angular momentum 
around the Sun-planet axis. We thus generalize here their work to a full 3D calculation, in 
order to investigate the influence of the radiation pressure on the trajectories (this paper), as 
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well as the density profiles and escape flux (following wo rks') 


This problem is similar to the so-called Stark effect ( Stark . 


19141 ): the effect of a constant 


electric held on the atomic Hydrogen’s electron. Its study can be transposed to celestial me¬ 
chanics in order to describe the orbits of artificial and natural satellites in the perturbed (e.g. 
by the radiation pressure forc e) Two -B ody Proble m. A recent description of the Stark effect 
solutions was already given by Lanto ine and Russelll (2 0111) . However, they give the analytical 
solutions of all trajectories only in the 2D case (and for bounded trajectories in the 3D case), 
and their formulas have so me issues as will be discussed later. Another analytical study pro¬ 
posed by IBiscani and Izzoi (20141) uses the Weierstrassian formulations to solve the motions for 


bounded and unbounded trajectories and to find periodic motions. Also, the motion can be 


approached numerically by develo ping the equations of motion in Taylor series 


jut this leads to 


some issues for high eccentricities (IPellegrini et al.l. 120141) . lHatten and Russelll (120141) compared 
recently these methods and their computing efficiencies. 


In this paper, based on the same formalism as 


Bishop and Chamberlain (119891) . we provide 


for the first time the complete exact 3D solutions of the Stark effect (and its celestial mechanics 
analogue) for any initial condition and for both bounded and unbounded trajectories. 

The first section describes the formalism used, before the sections [2j/[3j/|4] provide the equa¬ 
tions of motion and time. We then discuss about circular orbits in section [6], while a comparison 
with previous works is given in section 6, before we conclude in section [71 


2. Model 

In this work, we decide to study the effect of the radiation pressure on atomic Hydrogen in 
particular. Nevertheless, this formalism can be applied to any species subject to this force or to 
the interplanetary dust. We model th e rad iation pressure by a constant acceleration a coming 
from the Sun. According to Bishop! ( 1991 ). this acceleration depends on the line center solar 
Lyman-cc flux / 0 , in 10 11 photons.cm" 2 .s -1 .A -1 : 


a = 0.1774 /o (cm.s ) 

In spherical coordinates, the Hamiltonian of one Hydrogen atom can be written: 


( 1 ) 


+ P 2 e + v\ 

2m 2 mr 2 2 mr 2 sin 2 9 


( 2 ) 


GMm 


+ mar cos 9 
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with r the distance from the planet, 6 the solar angle, </> the angle with respect to the ecliptic 
plane, p ri po and p<p the conjugate momenta. —GMm/r represents the gravitational potential 
and mar cos 9 the potential energy from the radiation pressure acceleration a. An example of 
trajectory of a H atom subject to the radiation pressure is given in the figure [0 



Figure 1: Example of the trajectory of an atomic Hydrogen in the Earth exosphere. Left panel: view from the 
side. Right panel: view from the Sun. 


This problem is similar to the classical Stark effect f Starkl. 191411 : a constant electric field 


(here the radiation pressure) is applied to an electron (here an Hydrogen atom) attached to a 
proton (here the planet). Both systems are equivalent because the force applied by the proton 
(the planet) to the electron (the Hydrogen atom), i.e. the electrostatic force, varies in r -2 as the 
gr avitati onal force from the planet on the Hydrogen atom. Thus, we adopt the same formalism 
as iSo mm e rf eld] (119341 1 and use the parabolic coordinates. We use the transformation: 


u = r + x = r(l + cos6 l ) 
= r(l — cos 9) 


w 


r — x 


( 3 ) 
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with u and w always positive. Consequently, the Hamiltonian becomes: 


'H(u,w,p u ,p w ,p ( j ) ) 


2 up 2 . + 2 wp 2 p\ 2 GMm u — w 

___ _ _|-:_|_ uia _ 

m[u + tu) 2 muw u + w 2 


independent from t and 0. 


According to hamiltonian canonical relations, we have: 


Pu 


m(u + w) dw 
4 u dt 


Pw 


m(u + w) din 
4w dt 


(4) 


(5) 


d0 

Pit, = muw— 

^ dt 

2.1. Constants of the motion 

In this new system of coordinates, we study this Hamiltonian. First, EL is independent from 
t explicitly. EL is a conserved quantity along the time and corresponds to the mechanical energy 
E of the system. Moreover, as EL is independent from 0, according to canonical relations: 


dp^ _ d EL 
dt d0 


( 6 ) 


Thus, p^ is another constant of the motion. Once E and p<f, defined, the equation [4] can be 
rewritten: 

2 muE — 4 upi — — — m 2 au 2 + 2 GMm 2 
u 

(7) 

= —2 mwE + 4 wp 2 + — — m 2 aw 2 — 2 GMm 2 

w 

The left hand side is a function dependent only on u and p u , the right hand side depends 
only on w and p w . As both functions are equal and independent, they are equal to a constant 
A, a separation constant: 

A = 2 muE — 4 up 2 — — — m 2 au 2 + 2 GMm 2 

u 


= —2mwE + 4 wp 2 + — — m 2 aw 2 — 2 GMm 2 

w 


( 8 ) 
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The motion possesses three constants: E, A and p^. The equation 0 allows to express p v 
(respectively p w ) as a functions of E, A, p$ and u (respectively w): 


Pu = ± 


p w = ± 


-PM 

4 u 2 


QM) 

Aw 2 


(9) 


with 


P${u) = mau 3 — 2 mEu 2 — (2 GMm 2 — A)u + 


( 10 ) 


Qz{w) = maw 3 + 2 mEw 2 + (2 GMm 2 + A)w — p ^ 

2.2. Effective potentials 

We have already introduced the Hamiltonian EL of the system. We can extend the approach 
according to Hamilton-Jacobi equations: 

dS 


dqi 

dS 

dt 


= Pi 


= -EL 


( 11 ) 


where S is the Hamilton’s principal function or action. This function depends on initial con¬ 
ditions (as uo, wo, 0o and 1 0 ) and the actual position of the particle (as u, w, 0 and t). As 
previously demonstrated, EL = E and p$ are constants. Thus, 


as 

00 


p$ 


( 12 ) 


as = _ F 

dt 


and leads to: 


S = —E(t - t 0 ) +p</,(0 - 0o) +S[u 0 ,Wo,u,w,E,pf\ 
with S the part of the action independent from to, t, 0o and 0. 


(13) 


Moreover, the action S can be separated into two parts: one written as a function of u and 
p u coordinates, the other one with w and p w . By definition, according to the Hamilton-Jacobi 
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equations, we have: 


dS 

du 


as 

du 

as 


= Pu 


(14) 


= Pu 


as 

dt dw 

with p u (resp. p w ) a function only of u (resp. m), assuming E , A and p$ values already hxed 
by initial conditions. Then, we can separate again the action, leading to: 


S[u,w,u 0 ,w Q ,E,p^\ = S u [u,E,A,p,p] + <S™[m, E, A,p<f\ 

dS u 


(15) 


du 

dS w 

dw 


= Pu 


= Pu 


(16) 


According to the equation [9J. we have the following relations: 


d5u 

du 


E — 


pi GMm 


2 mu 2 


m 
~2 

m 

-{E-V u {v))> 0 


u 


A 

2 mu 


mau 


(17) 


d S w 
dm 


™ ( E Pi 

2 \ 2mm 2 

771 

-(E-V w (w))>0 


GMm A maw 
+ -+ x + —— 


m 


2mm 


with 


V u (u) 

pi 

GMm 

2 mu 2 

u 

V w (w) 

P\ 

GMm 

2mm 2 

w 


A mau 
2 mu 2 

A 


(18) 


maw 


2 mm 2 

V u and V w are effective potentials applied in u and m directions (represented in the figure [2]). 
These potentials play key roles for the motion because they constrained the motion in u and 
m directions independently. For the motion of the particle, we must respect two conditions: 
E > V u {u) and E > V w (w). These conditions are analogous to: 

P 3 (u) < 0 and Q%(w) > 0 (19) 

These both conditions are more restrictive than the usual E > E p where E p is the potential 


energy. 
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Figure 2: Representation in dimensionless unity of the shape of both potentials Vjj (corresponding to Vy, blue 
line) and Vw (corresponding to Vw, red line) for a set of E, A and p<f, values. The motion is possible only 
in the area where the potential is below the mechanical energy E, represented by the black horizontal line. 
The different roots of P 3 and Q 3 are displayed and correspond to the intersection of the potentials with the 
horizontal black line. U _ is the dimensionless value referring to U- (cf. section [2~T1 tabled]). Notice that Vw 
will cross only once the energy level £ for too high or too low £ values. 

2.3. Study of P 3 

P 3 is a polynom of degree 3 with lim P 3 (w) = + 00 . This polynom possesses three roots, 

W—>■+OO 

whose one is real at least. As P 3 (0) = > 0, one of these roots is real negative, according 

to intermediate value theorem since lim P 3 (w) = — 00 . Nevertheless, the motion occurs for 

u — y—OO 

positive u values, and we know this motion exists. It implies there is an interval in M + such 
as P 3 < 0 (otherwise, there is no motion in u-direction, no possible physically, cf. eq. fThlh To 
comply with this last condition, both other roots are real and positive. In summary, P 3 has 
three real roots: one negative and two positive. 

We call each root uq, U- and u + such as uq < 0 < w_ < u + and the w-motion is restricted 
to u G [«_;«+] (U G [[/_; U + \ in term of dimensionless quantities, as can be seen in the figure 
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2-4- Study of Q 3 

Q 3 is a polynom of degree 3 with lim Q${w) = + 00 . This polynom possesses three roots, 

w—t+oo 

whose one is real at least. As (£ 3 ( 0 ) = —p^ < 0 , one of these roots is real positive, according 
to intermediate value theorem since lim Qs(w) = + 00 . Nevertheless, the motion occurs for 

w —^-|-oo 

positive w values. We have restrictions on both other roots: they must be both real positive, 
both real negative or both complex conjugates. 

In the case where the three roots are real positive, we call each root wq, W- and w + such as 
0 < w_ < w + < w 0 and the motion is restricted to w G [w_; w + \ U [w 0 ; +00 [ (as can be seen in 
the figure [ 2 ] with the dimensionless quantities, cf. section [276]) . 

In the case with one positive root and both other complex or real negative, we call each 
root wq (the positive one), w_ and w + (keep the same order as previously defined if they are 
real) such as the motion is restricted tow 6 [w 0 ; +00 [ only ([Wo; +00 [ in term of dimensionless 
quantities). 

2.5. Restriction on the motion 

Each constant value of u or w defines a paraboloid in three dimensions. For each interval, 
constrained by fixed values of u and w, the motion will be contained between the paraboloids 
defined by these limit values as shown in the figure [3j 



Figure 3: Representation of the available space for the motion. Left panel: the blue area corresponds to the 
allowed region for the u-motion. Right panel: the red area corresponds to the allowed region for the w-motion. 
The regions are limited by paraboloids (parabolas here by projection). The Sun location is shown with a yellow 
circle. 

For the w-motion, this is always limited by two paraboloids defined by u = w_ and u = u + 
as shown by the blue area in the figure [3] (left panel). Similarly, for the w-motion, there are two 
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cases: the motion is constrained between one paraboloid (w = w o) and the infinity or between 
two paraboloids (w = w _ and w = w + ). Both cases are represented by the red area in the 
figure [3] (right panel). 

2.6. Summary of restrictions 

As previously demonstrated, the motion is constrained in specific areas of the 3D space. The 
motion is always between two paraboloids due to restrictions on u but it can be constrained 
between two other paraboloids or only one regarding w. Thus, the motion is constrained by 
four paraboloids or three, opened to infinity as shown on the figure [4] by the green area. 



Figure 4: Left panel: Combination of both panels in the figure [3] The final allowed region is in green. Right 
panel: representation in green of the trajectory in the figure [T| in the ( x, p ) frame at Earth. 


Nevertheless, there is at this point no other information on the exact motion of the particle. 
As shown in the figure [4] (right panel), the particle seems to explore all the area when its 
motion is restricted by four paraboloids. This is simply an observation. How can we prove 
that without any information on the exact motion of the particle? According to the Poincare 
recurrence theorem, if the dynamical trajectory of the autonomous system evolves in a finite 
volume of the phase space, then in any domain, as small as it could be, there are at least two 
points which belong to the same trajectory. Here, the motion is constrained in space with the 
four paraboloids but also in velocity because p u and p w are finite values and p$ is constant. All 
the positions in this part of the phase space can belong to the same trajectory. This is also 
linked to the KAM theory: here, the perturbation (the radiation pressure) affects the periodic 
motion (the ellipse, bounded trajectories) but it can remain quasi-periodic. As we will see, 










2 MODEL 


11 


the global motion is not periodic but u and w motions possess their own period according to 
another parameter. The global motion can be periodic only if all periods are commensurable. 

This is an important result because if the particle belongs to an area crossing the exobase 
and if this area is closed in the phase space, along a finite time, the particle will again cross the 
exobase. We can now define the ballistic and satellite particles as follows : both populations 
have no elliptic trajectories due to the radiation pressure, but they evolve in a closed domain. 
Depending on their constants of the motion, we can easily determine whether they cross (i.e. the 
domain crosses) the exobase or not, corresponding to ballistic and satellite particles respectively. 
Escaping particles are in the case where the initial value of w is higher than the highest real root 
of Q 3 and their available area is opened to the infinity. Thus, the theorem cannot be applied 
here. Even if the restriction area crosses the exobase, the particle may come from infinity, come 
close to the exobase and go away without crossing the exobase. We need, in order to identify 
escaping particles (that cross the exobase and go to the infinity), to track along the time the 
particle to know if these particles come from the exobase or not. Thus, it is necessary to solve 
their trajectory along the time. 


2.7. Dimensionless expressions 

For convenience, as usual in fluid mechanics, we define characteristic quantities. We decide 
to write all equations with dimensionless parameters. First, for distance, we define: 


n 


IGM 


■pressure 


( 20 ) 


This characteristic value was introduced by 


Bishop (119911 ) and defines the limit distance where 


the radiation pressure overwhelms the gravitation of the planet. Then, we rewrite: 


u 


UR 


■pressure 


( 21 ) 


lu vv i l pressure 

where U and W are the dimensionless quantities associated to u and w. The energy E is 
adimensionned with the use of the characteristic energy kBT exo . For p u and p w , we express 
them in units of \/mkBT exo , whereas p$ is expressed in units of \ZmkBT exo R P ressure- Finally, we 
choose mkBT exo yJGM / a as the unit for A and the unit for the adimensionless time r is derived 
from the units of u/w and p u /Pw We summarize these changes in table Q] 

We express all the previous equations as a function of these new quantities: 


£ 


2 UP?J + 21 VPb p .* 2A„ y _ 

U + W 2 UW U + W 2' 1 


( 22 ) 
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dimensional 

parameters 

unit 

dimensionless 

parameters 

u and w 

\J G Adj & R pressure 

U and W 

p u and p w 

\/ TYlk bR ex o 

P B and P w 

Pt 

\f mk B T exo GM/ a 

P* 

E 

^B^exo 

£ 

A 

mk B T exo ^/GM/a 

A 

V u and V w 

kB^exo 

Vu and Vw 

t 

1 GMm 

V ak B T exo 

T 


Table 1: Compilation of the transformations of the parameters into dimensionless ones 


with A a : 


p 2 

A = 2£U — AUPy — jj- + 2A a — X a U 2 

p 2 

= -2£W+ AWP 2 v + ^-2X a -X a W 2 

P 3 (U) = X a U 3 — 2£U 2 + (A — 2\ a )U + P 2 

= X a (U — U 0 )(U — UJ)(U — U + ) 

Qs(W) = X a W 3 + 2£W 2 + (A + 2X a )W-P 2 
= X a (W-W_)(W-W + )(W-W 0 ) 


Vu(U) = 

P 2 

2X a — A 

2 U 2 

2U 

V w (W) = 

p 2 

r <t> 

2X a + A 

2W 2 

2 W 


+ 


A a = 




am 


GMm 


knT 


kg T Rpressure 


X(Rpressure) 


(23) 


(24) 


(25) 


(26) 


the Jeans parameter at the distance R pre ssure■ 
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We can also reduce the equations of the motion. We introduce the dimensionless time r. 




'd U\ 2 

4 P 3 (U) 

.dr ) 

( u + w ) 

dW\ 2 

4 Q 3 (W) 

dr J 

(U + W) 2 

d 0 

P* 

dr 

uw 


(27) 


3. Dynamical trajectories 

In this part, we give implicit expressions for the dynamical trajectories of the particles, 
under the influence of both gravity and radiation pressure. Such expressio ns were a l ready g iven 


for t he 2D case with P^ = 0 and partially generalized to the 3 D cas e in 


(2011) with Jacobi elliptic functions (only bounded) (jJacobi, 


Lantoine and Russell 


1829) and in 


Biscani and Izzo 


(12014 1 with the Weierstrass functions for bounded and unbounded trajectories. This last work 
dealt with the dynamical trajectories of artificial satellites but it can be applied to exospheric 
species subject to the radiation pressure. We propose here corrections as well as a better way to 
give “simple” expressions for dynamical trajectories. We also provide expressions for unbounded 
trajectories that are missing in the literature. In the same way, we introduce our notations for 
the next papers to be published, where the influence of the radiation pressure on the density 
profiles and escape flux will be investigated. 

According to the previous part, we have different restrictions on the motion and thus, we 
must distinguish the cases, that will correspond to the different types of possible trajectories. 
We may thus define the ballistic/satcllite/escaping populations based on the roots of the P 3 
and Q 3 polynoms. The trajectories are not elliptic or hyperbolic at all but one can keep their 
basic definitions: crossing twice the exobase for ballistic particles, orbiting but not crossing the 
exobase for satellite particles, coming from the exobase and escaping to the infinity for escaping 
particles. 

The number of real roots is a key parameter for the analytical resolution of the trajectories. 
The previous equations of the motion, expressed as a function of time, must be rewritten as a 
function of a new variable % (the subscript 0 is necessary to distinguish from r) defined as: 


(U + W) dT 0 = dr 


(28) 
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Until we know the expression of U{%) and W{ 7o), we cannot express solutions as a function 
of time. The system of equations leads now to: 

2 


d U 

dTn 


dW 

d7o 


d (j) 

dTo 


-4 P 3 (U) 


VUW) 


PA h + w 


(29) 


dr = (U + W) d7o 
All solutions, trajectories and time, are parametrized according to the variable 7o without 


physical sense. 


3.1. The U-motion 

We solve the differential equation describing the u-motion: 

=-iPi(V) = -iK(V-U*)(V-UJ)(V-U + ) (30) 

with Uq, U- and U + being real roots. 

First, we set Y 2 = U — Uq. The motion occurs always with U > 0 and Uq < 0 to justify this 
change. We obtain: 

-AX a (U-U 0 )(U-U_)(U-U + ) 

4U 2 (|Q = -4A a Y 2 (Y 2 -U_ + U 0 ) (Y 2 -U + + U 0 ) 

<0 <0 

( 31 ) 

(W) = -UY 2 -U-+U 0 )(Y 2 ~U + + U 0 ) 

~{Y 2 -U. + U 0 )(Y 2 -U + + U a ) 
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Now, we set Z = 


V 


VU+-u 0 ‘ 

(-¥-\ 

\VKdrJ 


= (1-Z 2 )(Z 2 (U + -U 0 )-U. + U 0 ) 


(32) 


Finally, we have: 


dZ 


VKVU+ - UodTo 


= (1-Z 2 


\ 


U.-U 0 
Ua_ — Un 


<1 and >0/ 


We define k v as: 


1 U+-U- 1 ,.2 

u+ -U 0 U + - U 0 

<1 and >0 


U 


ku = 


U+ - U- 

u+-u n 


The final equation is: 


(33) 


(34) 


(35) 


dZ 


,^T a VU + -U 0 d% 

The solution of this equation is: 


) = (1 - Z 2 ) (Z 2 ~ (1 - kf,)) 


(36) 


z = MYKJuI - «,(T„ - Tv), k u ) (37) 

dn is a Jacobi elliptic function and Tu depends on initial conditions. 

The final expression for U is: 

' U(%) 

= U 0 + (U+ - U 0 )dn 2 (y/T a y/U + -U 0 (% - Tu), ku) 

< (38) 

= u+ - (U+ - U-)sn 2 (VKVU + -U 0 (T 0 - Tu), ku) 

^=U. + (U + - U_)cn 2 (VKVU+ - Uo(T - Tu), k v ) 

with cn and sn other Jacobi elliptic functions. These functions are 4iF(fc[/)-periodic (see Ap¬ 
pendix, eq. EQ). To determine Tu, we preferably choose the second expression for the sign. 
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The initial conditions give us 17(0), then: 


iarcsn 


! U+ - U(0) 
i u - U- 


, ku ) — \f~K l \/U + — UqTu 


(39) 


where arcsn is defined as the reciprocal function of sn: 

1 


arcsn(x, k) = 


dt = F(arcsin(x), k) 


lo Vi — t 2 V l — k 2 t 2 

with F the incomplete elliptic integral of the first kind (see Appendix, eq. IA.1D . 
The sign depends on Pu~. 


(40) 


d U 
dTo 


= 4 UPu 


(41) 


If Pu is positive (resp. negative), U increases (resp. decreases) and —sn 2 increases (resp. 
decrease s) too. Then, Tu i s pos itive (resp. negative). This solution corresponds to the solution 77 


given by 


Lantoine and Russelll (j2Qll|). The [/-motion is always constrained. The characteristics 


of the motion/particle (ballistic, satellite or escaping) is thus determined by conditions on 

Qs(w). 

3.2. The W-motion 

For the w-motion, we need to solve the differential equation: 


d W 


= 4 Q Z {W) = JX a (W - W 0 )(W - WJ){W - W + 


(42) 


Wo is a positive real root but we must distinguish between various cases for the other roots: 
W+ and W- may be real or complex conjugate roots. 


3.2.1. Three real roots 

We must consider two cases: W > W 0 and W + > W > W-. The second one occurs 
mathematically and not only physically when and W + are positive. 

a) For W_ < W < W + . This case occurs for bounded trajectories (ballistic or satellite trajec¬ 
tory). Here, this is the same treatment as for the [/-motion. After setting the two following 

Y 

substitutions Y 2 = Wq — W and Z = —=====. we obtain a similar solution: 

y/Wo — W- 


Z = dn(VVVWo - Win - Tw ), k w ) 


(43) 
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The final expression for W is 

' W{%) 


= W 0 - (bF 0 - - Tw), k w ) 

= W+- (W + - W.)cn 2 (^Y a ^Wo-W4To - Tw), k w ) 

^ = W. + (W+ - W.)sn 2 (y/T a VW 0 - W-(% - Tw), kw) 
The initial conditions giving us W (0) and then leads to: 


(44) 


iarcsn 


The sign depends on P w \ 


= 'NJw^W-Tu 


d W 

d% 


= 4WP] 


w 


(45) 


(46) 


If Pw is positive (resp. negative), W increases (resp. decreases) and sn 2 increases (resp. 
decreases) too. Then. Tw is n egat i ve (re sp. positive). This solution corresponds to the solution 


£/ given by 


Lantoine and Russell (2 01 if ) 


b) For W > Wo . First, the motion occurs only for W E [lYo; +oo[ (corresponding to escaping 
particles). Then, we set Y 2 = W — W 0 - We obtain: 


dW 
d To 


4A a (W - W 0 )(W - W_)(W - W+ 


4 Y 2 


= 4A a Y 2 (Y 2 -W_ + Wp) (Y 2 -W + + Wp) 
>o >o 


dY 

dTo 


\ a (Y 2 -W- + W 0 )(Y 2 -W + + Wo) 


(47) 


dY y 




(Y 2 -W- + W 0 )(Y 2 -W+ + Wo) 


Now, we set Z = 


Y 


VW 0 - W- 
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Finally, we have: 


dZ 


v'WWo - W_dT 0 


= (Z 2 + 1 ) 


2 , W 0 - W 4 




W 0 - W- 

<1 and >0 


We define kw as: 


w 0 -w + _ w+-w- 2 

W n - W- W n - W_ w 


kw = 


<1 and >0 


' w + - w_ 

Wo - W_ 


The solution of this equation is 


Z = cs(ytvw„- W_(To - Tw), kw) 


(48) 


(49) 


(50) 


(51) 


cs is a Jacobi elliptic function and Tw depends on the initial conditions. 

The final expression for W is: 

' W(To) 

= Wo + (Wo - W.)cs 2 (VKVWo- W_(T 0 - Tw), k w ) 

< (52) 

= W + + (W 0 - W_)ds 2 ( v W:VH / o - W_(T 0 - T w ), k w ) 

k = W_ + (Wo - W-)ns 2 (VKVWo-W4% - T w ), k w ) 

These functions are 4A"(/qy)-periodic and are defined on % — Tw £ M/{4miF(fc w )|m 6 Z}. 
These functions diverge at 4mK(kw) but this is not an issue: the motion can diverge with 
respect to %, but, according to the time r and the integration (see Section HD, the re-motion 
remains continuous for r G R as 


To — Tw e]0; 4A (kw) t G M. 


Let us assume the initial conditions provide W(0), then: 


f / W(0)~ Wp 
\\ W 0 -W_ 



vW^o - W.Tw 


iarccs 


(53) 
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with 


arccs(x, k) = 
The sign depends on P w : 


f , , = dt — F ( arctan I — ) , k 

I o VT+CVt 2 - 1 + k 2 V U ' 


d W 
d7o 


= mp w 


(54) 


(55) 


If Pw is positive (resp. negative), W increases (resp. decreases) and cs 2 increases (resp. 
decreases) too. Then, Tw is positive (res p. negative). This solution corresponds to the solution 
£// given bvl Lantoine and Russell (2 0111 1. These expressions are useful only for three real roots. 
A last case remains: only one real positive root. 

3.2.2. Only one real positive root 

For any initial conditions, the particle will be escaping if Q 3 has only one real root. The 
solution for thi s case is more complex compared with the previous solutions. As done by 


Lantoine and Russelll (1201 11 1, we could apply some transformations and obtain a new expression, 


that would be a combination of cn and sn. Nevertheless, we propose a simpler way to determine 
the expression with the knowledge of all roots, even imaginaries. We start from the equation 
[561 with W + and W_ not real: 

2 


d W 
d% 


= 4A a (W - W 0 )(W - W+){W - W_) 


(56) 


After the separation of variables, we obtain: 

f w _dlT_ 

J \J(W' - W Q ){W' - W + )(W' - WJ) 

Now, we apply the procedure proposed in 


— 2v / Ai(7o — Tv 


w) 


(57) 


Abramowitz and Stegunl (119641) (p. 597) in the 


case where we have only one real root. First, we define: 

A 2 = \J (IT 0 - W+)(Wq - W-) = V(7F 0 - Re(IF +)) 2 + Im(IF + ) 2 = 
and also 


^3(^0 

A a 


(58) 


kw 0 = A I 7 


1 1 W 0 - W+ + W 0 - W_ 


A 2 


1 1 Wo - Re(W + ) 

2 ~ 2 ~ 


A 2 


(59) 


According to 


Abramowitz and Stegunl (1196411 (p. 597), the left hand side corresponds to: 

f w dW 1 F(9, k 


^Wo) 


I Wo V( w ' - w o){W’ - W+)(W - W-) 


A 


(60) 
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with F the Elliptic function of first kind (defined in the appendix Appendix A ) and 

A 2 — (W — W 0 ) .. , 0 1 — cos 9 


C0S ° A 2 + (W-W 0 ) 

According to the equation 1571 


W = Wq + A 2 


1 + cos 6 


(61) 


6 — am(2A\/A^(7o — Tw 0 )i kw 0 ) 


(62) 


with am the Jacobi amplitude defined as 


am(x, k) = F 1 (x, k ) 


(63) 


Thus 


W(T 0 ) = W 0 + A 2 


1 — cn(2A\/A^(7o — Tw 0 ) , kw 0 ) 
1 + cn(2Av / ^(7o — Tw^ r kw rr 


Finally, one can simplify the expression based on 


Olver et al. 


(64) 


( 20101 ) (equation 22.6.18): 


W{%) = W 0 + A 


2 dn 2 (A-\/A^(7o — Two), kw 0 ) 


(65) 


cs 2 (Aa/A^(7o - Tw 0 ) , kw 0 ) 

This function is IK (fcy^ 0 )-periodic and is defined on %—Fw £ {2K{kw 0 )-\-IrnK(kw 0 )\rn E 

Z}. This function diverges at 2 K(kw 0 ) + 4 mK(kw 0 ) but this does not impact our result: the 
motion can diverge with respect to 7o, but, according to the time r and the integration (see 
Section HD, the re-motion remains continuous for r 6 1 as 


% — Tw 0 e] — 2A (k w0 ); 2 K(k Wo )[<^=^- r E 
As usual, according to the initial conditions, we define 7 w 0 as 

'A 2 - (W(0) - W 0 ) 


with 


iarccn 


arccn(x, k) = 


A 2 + (W(0) - W 0 ) 


, k Wo I = 2Av / Ai'7v 


Wo 


df = F(arccos(x), k ) 


( 66 ) 


(67) 


y/l + t 2 \/1 — k 2 + k 2 t 2 

If Pw is positive (resp. negative) then 7ux is negative (resp . positive). This solution 
corresponds to the solution 77 given by Lantoine and Russell (12011 1 but our derived solution is 
less complex to use (with less time computing). 
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4. Time equation 


We gave analytical formulations for the different kinds of trajectories, expressed implicitly. 
Now, since we have all expressions of the trajectories as a function of To, we can express the 
real time r (or t ): 


t = 


-r 


GMm 
ak B T e , 

t(%) = [ V(T') + W{T') dT' 


( 68 ) 


Here, we use MATHEMATICA and MAPLE to derive the primitives. It is necessary to be 
very careful since these different programs can have different definitions of the elliptic functions 
for example. To avoid these issues, we remind at each use the definition employed. The first 
part of the integral gives: 

[•To 

/ U(T')c\T' 


= U 0 % 


(69) 


(X 

+— [E(am(a(% - %), k n ), k v ) - E(am(-aTu, %),%)] 

^ a 

with a = y/X^y/U + — U 0 , E the incomplete elliptic function of the second kind (see Appendix, 
eq. IA.II) . 

The second part of the integral is more complex because we have different expressions 
according to the number of roots and the initial conditions. In the case of three real roots, if 
W+ >W( 0) > then 

[•To 

/ W{T)dT 
Jo 


= W 0 T 0 


(70) 


(3 


A 


a 


[E( am(/3(7o — 7w), k w ), kw) — E(am(—f3Tw , kw), kw} 


with /3 = \JXJVWo - W-. 
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If W (0) > Wo with three real roots then 



W{T) dT' 


= W 0 To 


P 


(71) 


[E(am(P(T 0 - Tw), kw), kw) ~ E(am(-/3Tw, k w ), k w )] 


P 


en(/3(7o — Tw), kw) cn(— (3Tw, kw) 


A a lsd((3(% — Tw), kw) sd(— PTw, k w ) 

Finally, in the case of only one real root, the time equation is given by: 

[•To 

/ W(T)dT 


= (Wo + A 2 )To 


”7 

— — [if(a m ( 27 ( 7 o — Tw 0 ), k Wo ), k Wo ) — E(am(—2^Tw 0 , k Wo ), k Wo )] 


+ J_ sn(27(7o - Two), fcw 0 )dn(27(7o - T Wo ), k Wo ) 

K [ 1 + cn( 27 ( 7 o — 7w 0 ), kw 0 ) 


sn (-2-fTwo, k Wo )dn(—2'y7w 0 , k Wo ) 


1 + cn(— 277 vvo> k Wo ) 


72) 


with 7 = Aa/A^- 

For both last cases, % — Tw belongs to ]0]AK(kw)[ (respectively ] — 2K(kw 0 )]2K(kw 0 )[) 
for the equation 1711 fresp. for the equation [72]) . Nevertheless, when To — Tw tends to 4 K(kw 0 ) 
(resp. 2K{kw 0 ))i the integration diverges and r too. The w- motion occurs on an subset of 
M with respect to To but on M entirely with respect to r. The transformation of r into 7o is 
bijective because the integrand U + W is strictly positive. 


5. ^-equation 


To complete the description of the motion as a function of time, it is also necessary to solve 
the evolution of the angle 0, obeying to: 

d0 


dT 0 


P * \U(%) + l V(%)) 


(73) 


r-To 


0(T O )-0(O) = / P[ 


1 1 
+ 


U(V) W(V) 


d T 
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As already done in the previous part, we separate into two integrals. The first part still gives: 

rT 0 x 

Jo Wr) 


d T 


n ( 1 — ain(ft(To - Tu),k v ), k v 
aU+ 


(74) 


n ( 1 — ^j--,am(-aTu,k n ),k„ 
all. 


with II the incomplete elliptic integral of the third kind (see Appendix, eq. IA.1D . 


In the three real roots case for Q 3 , if initially we have W + > hh(0) > HA: 

f V 1 

Jo W^) dT0 

n ^1 - am (P(T - Tw),k w ), kw^j 

= pwl 


n 


w+ 

HA 


; am(— PTw, kw),k 



PW- 


or if initially W( 0 ) > W 0 : 



d T 


P%-u 


HA 

W 0 - HA 



PW- 


(75) 


(76) 


n 


HA 

W 0 - HA 


; am(~PTw, kw ), k 



pw_ 
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and for the last case with one real root: 

r To 1 

Jo Wr r ) 


dV 


To 


o 


Wo - A 2 


(W 0 + A 2 


n 


4W 0 7(W 0 - A 2 ) V 4A 2 W C 


(Wq - A 2 ^ 2 


am( 27 ( 7 o — 7w 0 ), kw 0 ), kw 0 


i 1 arrtan ( p <t> sn ( 27 ( 7 o — Tw 0 ), k Wo ) 
2 \ZyWoMMTo-Two),kwo) 


(77) 


(Wo + A 2 ) (Wo - A 2 ) 2 , 0 ^ , , 

-n I- 7~ttt7~7z —; am(—2^Tw 0 , k Wo ), k Wo 


4Wq7(W 0 - A 2 )' 


4 A 2 W 0 


1 ( Ptj, sn(— 277 wq, fcyy 0 ) \ 


2Prb V 2 7 W 0 dn(—2 ^Tw 0 i k Wo ) J 


This last formula is only available for 7o — 7iv 0 s] — 2K(k\y 0 ); 2K(kw 0 )[, range where it is 
continuous. 


6. Circular orbits 

With the solutions previously derived, we can know the exact motion of a bounded or un¬ 
bounded particle as a function of the time such as given in figure [6j It is clear that even a 
bounded trajectory the motion has no periodicity at all (especially for the 0 motion). Nev¬ 
ertheless, it coul d be in teres tin g to focus on stable bounded orbits and search for periodic 
motions (as in Biscani and Izzol (1201411 ). and thus investigate in particular the circular stable 


orbits for spacecrafts flNamouni and Guzzol. 


2007) or the possible positions for satellite particles 


produced by collisions in the exosphere flBeth et all 2 0141) . Thus, we dedicate this section to 
the conditions to obtain such orbits. 

For a specific set of initial conditions, it can be possible to obtain circular orbits. This orbit 
occurs when, on the one hand, the attraction of the planet projected along the x-axis is equal 
to acceleration due to the radiation pressure: 

GMmx . 

--- ma = 0 (78) 


— ma = 0 
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In dimensionless quantity, this will be expressed by: 


R 2 + cos 0 = 0 


(79) 


On the other hand, it is also necessary that the centrifugal force induced by the rotation 
around the x-axis is equal to the acceleration around the planet in the perpendicular plane to 
the x-axis. Thus, we obtain the secondary equality: 


GMmp + mv% ^ 


P 


In dimensionless quantity, this will be expressed by: 


\ a Rsm 4 9 — Pi = 0 


(80) 


(81) 


Combining these two equations, we obtain: 


Pi 


(sin 2 0) 9 — (sin 2 0 ) s + —A = 0 

A n 

We need to study the polynom 


(82) 


P(X) = X s - X s + 


p 2 

K 


(83) 


with X = sin 2 0 6 [0; 1]. Depending on the P 2 values, we have zero, one or two solutions for 
P(X) = 0. Indeed, P(0) = P(l) > 0 so that according to the Rolle theorem, there is a e]0; 1[ 
with P\a ) = 0. Here, a is equal to 8/9. Nevertheless, if P(a) is positive then P(X) = 0 does 
not have solutions and inversely, if P(a) is negative then we have two solutions. This value is: 


p (-1 = !jl - 1 (- 
9 9 V 9 


A critical maximum value of P (l , thus exists to allow for circular orbits and is 

8 \/A^ 


(84) 


Ia*i = 


9^3 


(85) 


Above this value, we cannot find any bounded trajectories: there is no any equilibrium point 


and thus no circular orbits (stable or not). For lo wer val ues of |PJ , we have 


;wo solutions for 


P{X) = 0: one stable and one unstable as shown by Namouni and GuzzoJ fi2007h and theirs plots 


of the equipotentials. These two solutions correspond respectively to the stable point around 
which the equipotentials are closed and to the saddle point, which is the last limit where one 
can find closed equipotentials, and is the only point where two equipotentials can cross. As long 
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as IP^I < \P c <f>\, these both specihc points exist: they have the same U. Physically, the potential 
Vw has two extrema as plotted in figure [2] When \P^\ reaches \P c <t>\, the local minimum goes 
to the right and the local maximum goes to the left at the same location W cr u- For higher \P^\ 
values, Vw have no extremum anymore: the potential is strictly decreasing and the particles 
are unbounded (escaping). Thus, the bounded particles, satellite and ballistic particles, have 
|P<0 < |P C 0|. The distinction between them thus depend on if they cross or not the exobase. 

The critical value s are given in (z, p) coordinates (z is — x for us, in comparison with 
Namouni and Guzzol ( 20071 )). In dimensionless unities and in ( R,9 ) using the equations I8T1 
then 1791 the critical orbit is: 

2V2' 


(Rent, Ocrit ) = | VT - arCSlll 


or in (U, W ) coordinates: 


( U cr it, bberit) 


(- 


4 


( 86 ) 


(87) 


V3a/3’ 3^3, 

The real positive roots of the polynomial 1831 combined with the equality [79] give the positions 
of the circular orbits (two coordinates are necessary) allowed to spacecraft or particles under 
the influence of both gravity and radiation pressure. 


7. Summary 

The knowledge of the exact trajectories of particles or satellites under the influence of gravity 
and radiation pressure needs the calculation of the spatial coordinates, i.e. the U/W/(f> motions, 
as well as the time evolution. We summarize all needed equations in table [7] 

The U -motion is provided by the equation [38] The W -motion is provided by the equations 
1441 [52] or [65] The 0-motion is provided by P^ x [741+ (|751 or 1761 or ITT])]. The time equation is 
provided by [69]+ m or ED or 1721) ]. All the expressions are functions of 7o that is not the real 
time. We thus have implicit expressions as a function of time. The function r(7o) is bijective 
but cannot be inversed analytically, a numerical inversion is needed to derive the real time. 

Besides, our 3D solutions can be easily applied to the 2D case. Indeed, in the 2D case, P^ = 0 
and thus, one of the roots for each polynomial P3 and Q 3 is null: it could be Do or P_ for P3 
(if U + = 0, there is no possible motion) and any roots of Q 3 . We precise that in this case the 


0-motion is not important because the motion is planar. Compared with 


Lantoine and Russell 


((2011), our formulations are first developed for the 3D case and can be used easily for the 2D 
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cases 

Three real roots 

for P 3 (always) 


Three real roots 

for Q 3 

One real root 

for Q 3 

trajectory 

bounded 

bounded 

unbounded 

unbounded 

[/-motion 

IR-motion 

E1L2D) 

EKL2D) 

[52KL2D) 

Els) 

time equation 

E1L2D) 

EDDIED) 

EIFnew* 

E2Fnew* 

0 -motion 

[74KL3D) 

E5(L3D) 

ESFnew* 

EZFnew* 

7q range 

% E M 

% e M 

% G [0; (1 — sg(7w))2A (k w ) + Tw[ 

7q € [0; 2K(k Wo ) + ' 


Figure 5: Summary of the different solutions for each kind of motion. The notation *NEW* corresponds to 
the new solutions derived in this paper, not yet derived in previous works. The symb ol ( a) correspo nds to 
an example of formula with a simpler ex pression than the one propo sed by 


labeled solutions were explicitely given by 


solutions were also given by 


Lantoine and Russell! ( 201 11. L2D 


Lantoincan^Russcll ( 2011 1 only in the 2D case, whereas L3D labeled 


Lantoine and Russell (120111) in the 3D case. The function sign is noted sg. 



Figure 6: Plots of a bounded particle (or spacecraft) motion in the ( X , p) plane (upper left panel), of the time as 
a function of 7o (upper right panel), of the motion in polar coordinates (0, p) (lower left panel) and the U — W 
coordinates as a function of the time (red for U, blue for W, lower right panel). U and W do not show any 
periodicity because their periods are not commensurable (i.e. the ratio is a rational number) with 7o and thus 
with the time r 
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case, whereas 


Lantoine and Russell (2011J) gave only the methodology to obtain the 3D solu¬ 


tio ns based o n 2D ones but not the expressions, which apparently leads to complex expressions. 


Biscani and Izzo ('2 0141 ) provided also the exact formulas for boimded and unbounded trajec¬ 


tories using the Weierstrass functions but this formulation is also difficult particularly because 
of the need to use the Inverse Weierstrass function, not implemented in all computer softwares 
and the need to work with complex values (e.g. the complex logarithm function). In this paper, 
we solved the motion for the bounded and unbounded trajectories in the 3D case; we provide 
the exact formulas for all cases, as well as the definition s used in |Appendix A We al so highlight 
in table [7] which solutions are simpler compared with lLantoine and Russelll (120111). w hich are 


completely new and which were only provided in the 2D case by 


Lantoine and R usse ll (2 0111) . 


Moreover, beyond the new exact solutions given in this paper, the derivation of our sol utions 


based on Jacobi elliptic functions allows a good computing time and accuracy. 


Hatten and Russell 


compared thr e e type s of solutions for the Stark effe ct: two exact ones, pro posed by 


Lantoine and Russell (2011) (Jacobi elliptic functi ons) and 


strass elliptic functions), and a numerical one by 


Biscani and Izzol (120141) (Weier- 


Pellegrini et al 


( 20141 ) (based on Taylor se¬ 


ries) . They compared the CPU time, the numbe r of c alls for each analytic ellipti c function 


and the accuracy between 


Biscani and Izzol (120141) and lLantoine and Russelll (120111) . Even if 


we do not agree wit h t 


re number of evaluations of each Jacobi elliptic function mentioned 


by Hatten arid Russell ( 20141 ) (e.g. to call cn and sn is similar to call am then cos and sin: 


cn = coso am and sn = sino am), two arguments show our solutions are efficient in terms 
of CPU time and accuracy : first, t he soluti ons expresse d in terms of Jacobi elliptic functions 


(such as in this paper or by 
elliptic functions (used by 


Lantoine and Russell (2011)) are more efficient than Weierstrass 
- oanu - JiL- 

Biscani and Iz zol (I2014D) : seco n d, se veral solutions given above are 


less complex to implement than those by 


Lantoine and Russelll (2 01 lh . e.g. equation 1651 Also, 


the analytical formulations are preferable for long duration motions. 


8. Conclusions 

We determined analytically the trajectories of the particles or spacecraft under the influence 
of both planetary gravity and stellar radiation pressure. We thus provide for the first time the 
complete exact solutions of the well-known Stark effect (effect of a constant electric field on the 
atomic Hydrogen’s electron) with Jacobi elliptic functions, for both bounded and unbounded 
orbits. These expressions may be implemented for modeling spacecraft or particles trajectories: 
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instead of solving the equation of the motion, based on differential equations, with numerical 
methods such as the Runge-Kutta method where one cumulates errors along the time, it is 
here possible to obtain precise expressions of the motion with only periodic errors, due to 
the precision on the evaluation of the elliptic functions used. In particular, we provide the 
analytical conditions for stable circular orbits. Moreover, we discuss about the possible issues 
inherent to the formalism used and the importance of being extremely careful with the routines 
implemented. 


re formalism used here will allow us in a next paper to generalize the work by Bishop and Chamberlain 


(119891) to derive the exact neutral densities and escape flux in planetary exospheres, under the 


influence of both gravity and stellar radiation pressure. This is important for understanding 
the atmospheric structure and escape of planets in the inner solar system, as well as the at¬ 
mospheric erosion during the early ages where the radiation pressure (and UV flux) of the Sun 
was extreme. 
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Appendix A. Elliptic integrals 


In this paper, we use the three incomplete elliptic integrals F, E and II: 

1 


Ffak) = 


'o \/l — k 2 sin 2 9 


dd 


E(<t>, k) = / \/l-k 2 sm 2 6d6 


(A.l) 


11(77,; 0, k) = 


d0 


'0 (I — 77 sin 2 6) \Jl — k 2 sin 2 6 
Sometimes, other formulas (shown below) are proposed with the change sind = t but one 
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needs to be very careful: this change is bijective only for 6 G [—7 t/2 ; 7 t/2 ] 


““Vi - w 


^ cos 9 


i ,-- d t — . . 

/ 0 \/l - t 2 Jo I cos 6> I 


\/l — k 2 sin 2 6 d 6 


r'X=sin(f) 


: d t = 


^ cos 9 


Vl - k 2 t 2 V 1 - t 2 Jo I cos 6>I yi _ k 2 sin 2 9 


de 


(A.2) 


: dt = 


cos 9 


(1 - nt 2 )y/l - k 2 t 2 y/l - t 2 Jo I cos 6»| (i — n sin 2 9)\/l ^ k 2 sin 2 6 


d6 


These expressions are not exac tly E, F an d K 


the range (j) G] — 7 t/2 ; 7 t/2 [. L ant oine an d Russell 1201 ill did not precise which formulations they 


'hey 


agree with the previous formulas IA.1I 


m 


used. According to theirs formulas and results, they used the left-hand side of the equations 
IA.21 This may be a problem for bounded trajectories: for (ft = am(7o) x = sn(7o), the 
integrals IA.2I are not continuous contrary to IA.1I Depending on the computer software, the 
routines and the definitions used for these functions, the results can show some issues (e.g. no 
continuous motion). 
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